/*
 * 简单矩阵实现: 线性运算(+,-,数乘,数除), 初等行变换, 矩阵乘法, 转置, 求逆
 * FMatrix: 固定大小矩阵
 *
 * 历史:
 *	2019-12-15 姚彧 创建
 *  2019-12-16 ...  加入名字空间yy::math::
 *                  使用vs2019进行文件格式化
 *  2019-12-30 ...  名字空间改为qk
 */

#ifndef _QK_FMATRIX_HPP_
#define _QK_FMATRIX_HPP_

#define FMATRIX_TEMPLATE template <int m, int n, class Real, int size = m *n>
#define FMATRIX_T_TEMPLATE template <int m, int n, class Real, class T, int size = m *n>
#define FMATRIX_MN FMatrix<m, n, Real>

#include <cstring>
#include <cmath>
#include <algorithm>
#include <initializer_list>

namespace qk
{
	namespace math
	{

		template <int _m, int _n, class Real, int _size = _m * _n>
		struct FMatrix
		{
			typedef Real		value_type;

			Real v[_m][_n];
			int m;
			int n;
			int size;

			FMatrix() : m(_m), n(_n), size(_size) { makeZero(); }
			FMatrix(int _m0, int _n0) : m(_m), n(_n), size(_size) { makeZero(); }
			FMatrix(int _m0, int _n0, std::initializer_list<Real> values)
				: m(_m), n(_n), size(_size)
			{
				int i = 0;
				makeZero();
				for (auto iter = values.begin(); i < size && iter != values.end(); ++i, ++iter)
					v[0][i] = *iter;
			}

			// 二维下标运算符
			Real& operator()(int i, int j) { return v[i][j]; }
			const Real& operator()(int i, int j) const { return v[i][j]; }

			// 一维下标运算符
			Real& operator[](int i) { return v[0][i]; }
			const Real& operator[](int i) const { return v[0][i]; }

			FMatrix& makeZero()
			{
				memset(v[0], 0, sizeof(v));
				return *this;
			}

			FMatrix& makeIdentity()
			{
				makeZero();
				for (int i = 0, s = std::min(m, n); i < s; ++i)
					v[i][i] = (Real)1;
				return *this;
			}
		};

		//初等行变换
		FMATRIX_TEMPLATE
			void Eij(FMATRIX_MN& A, int i, int j)
		{
			for (int k = 0; k < n; ++k)
				std::swap(A(i, k), A(j, k));
		}

		FMATRIX_T_TEMPLATE
			void Ei_c(FMATRIX_MN& A, int i, T c)
		{
			Real c0 = c;
			for (int k = 0; k < n; ++k)
				A(i, k) *= c0;
		}

		FMATRIX_T_TEMPLATE
			void Eij_c(FMATRIX_MN& A, int i, int j, T c)
		{
			Real c0 = c;
			for (int k = 0; k < n; ++k)
				A(j, k) += A(i, k) * c0;
		}

		// 一元运算
		FMATRIX_TEMPLATE
			FMATRIX_MN operator-(FMATRIX_MN A)
		{
			for (int i = 0; i < size; ++i)
				A[i] = -A[i];
			return A;
		}

		// 线性运算
		FMATRIX_TEMPLATE
			FMATRIX_MN operator+(FMATRIX_MN A, FMATRIX_MN const& B)
		{
			return A += B;
		}
		FMATRIX_TEMPLATE
			FMATRIX_MN& operator+=(FMATRIX_MN& A, FMATRIX_MN const& B)
		{
			for (int i = 0; i < size; ++i)
				A[i] += B[i];
			return A;
		}

		FMATRIX_TEMPLATE
			FMATRIX_MN operator-(FMATRIX_MN A, FMATRIX_MN const& B)
		{
			return A -= B;
		}
		FMATRIX_TEMPLATE
			FMATRIX_MN& operator-=(FMATRIX_MN& A, FMATRIX_MN const& B)
		{
			for (int i = 0; i < size; ++i)
				A[i] -= B[i];
			return A;
		}

		FMATRIX_TEMPLATE
			FMATRIX_MN operator*(FMATRIX_MN A, Real k)
		{
			return A *= k;
		}
		FMATRIX_TEMPLATE
			FMATRIX_MN operator*(Real k, FMATRIX_MN A)
		{
			return A *= k;
		}
		FMATRIX_T_TEMPLATE
			FMATRIX_MN& operator*=(FMATRIX_MN& A, T k)
		{
			Real k0 = k;
			for (int i = 0; i < size; ++i)
				A[i] *= k0;
			return A;
		}

		FMATRIX_T_TEMPLATE
			FMATRIX_MN operator/(FMATRIX_MN A, T k)
		{
			return A /= k;
		}
		FMATRIX_T_TEMPLATE
			FMATRIX_MN& operator/=(FMATRIX_MN& A, T k)
		{
			Real k0 = k;
			if (k0 != (Real)0)
				A *= (Real)1 / k0;
			else
				A.makeZero();
			return A;
		}

		//转置
		FMATRIX_TEMPLATE
			FMatrix<n, m, Real> Transpose(FMATRIX_MN const& A)
		{
			FMatrix<n, m, Real> B;
			for (int i = 0; i < m; ++i)
				for (int j = 0; j < n; ++j)
					B(j, i) = A(i, j);
			return B;
		}

		//矩阵乘法
		template <int m, int n, int s, class Real>
		FMATRIX_MN operator*(FMatrix<m, s, Real> const& A, FMatrix<s, n, Real> const& B)
		{
			FMATRIX_MN C;
			for (int i = 0; i < m; ++i)
				for (int j = 0; j < n; ++j)
					for (int k = 0; k < s; ++k)
						C(i, j) += A(i, k) * B(k, j);
			return C;
		}

		//矩阵求逆
		template <int n, class Real>
		FMatrix<n, n, Real> Inverse(FMatrix<n, n, Real> A)
		{
			const Real one = 1;
			const Real zero = 0;
			int i, r, ri;
			Real Aii, c;
			Real determinant = 1;
			FMatrix<n, n, Real> B;
			B.makeIdentity();
			// 向下消元
			for (i = 0; i < n - 1; ++i)
			{
				//第(i,i)从第i列选主元
				Aii = abs(A(i, i));
				ri = i;
				for (r = i + 1; r < n; ++r)
				{
					if (Aii < abs(A(r, i)))
					{
						ri = r;
						Aii = abs(A(r, i));
					}
				}
				//换行元
				if (ri != i)
				{
					Eij(A, i, ri);
					Eij(B, i, ri);
				}
				//主元归一
				Aii = A(i, i);
				if (Aii != one)
				{
					determinant *= Aii;
					c = one / Aii;
					Ei_c(A, i, c);
					Ei_c(B, i, c);
				}
				//向下消元
				for (r = i + 1; r < n; ++r)
				{
					c = -A(r, i);
					Eij_c(A, i, r, c);
					Eij_c(B, i, r, c);
				}
			}
			//最后一行
			determinant *= A(i, i);
			Ei_c(B, i, one / A(i, i));

			//回代
			for (i = n - 1; i > 0; --i)
				for (r = i - 1; r >= 0; --r)
					Eij_c(B, i, r, -A(r, i));

			return B;
		}


		// 求行列式
		template<int n, class Real>
		Real Determinant(FMatrix<n, n, Real>A)
		{
			const Real one = 1;
			const Real zero = 0;
			int i, r, ri;
			Real Aii, c;
			Real determinant = 1;
			// 向下消元
			for (i = 0; i < n - 1; ++i)
			{
				//第(i,i)从第i列选主元
				Aii = abs(A(i, i));
				ri = i;
				for (r = i + 1; r < n; ++r)
				{
					if (Aii < abs(A(r, i)))
					{
						ri = r;
						Aii = abs(A(r, i));
					}
				}
				//换行元
				if (ri != i)
				{
					Eij(A, i, ri);
				}
				//主元归一
				Aii = A(i, i);
				if (Aii != one)
				{
					determinant *= Aii;
					c = one / Aii;
					Ei_c(A, i, c);
				}
				//向下消元
				for (r = i + 1; r < n; ++r)
				{
					c = -A(r, i);
					Eij_c(A, i, r, c);
				}
			}
			//最后一行
			determinant *= A(i, i);

			return determinant;
		}
	} // namespace math
} // namespace qk
#endif //_QK_FMATRIX_HPP_
